/*
 * td4.cpp
 *
 *  Created on: 13 févr. 2012
 *      Author: aelise
 */
#include "td2.hpp"
#include "td3.hpp"
#include "td4.hpp"
#include "bmp.hpp"

void writeFirstQuarter(double * final, double * quarter, unsigned int x, unsigned int y,  int py){

	int offset = 0;
	int offsetQ = 0;
	for(unsigned int i = 0; i < y; i++){

		memcpy(final + offset, quarter + offsetQ, sizeof(double) * x );
		offset += py;
		offsetQ += x;

	}

}

/*
 * Copie le premier quart d'une image dans un nouveau tableau
 *
 * m : image de départ
 * x : taille de m en x
 * y : taille de m en y
 *
 * retour : image 4x plus petite contenant le quart superieur gauche de m
 */
double* extractFirstQuarter(double* m, unsigned int x, unsigned int y){

	int offset = 0;

	int size = static_cast<int>(y*x/4);

	double* tmp = new double[size];

	for(unsigned int u=0; u< x*y/2; u += x) {

		memcpy(tmp + offset, m + u, sizeof(double)*x/2);
		offset += x/2;
	}

	return tmp;
}


//tmp = addFirstQuarter(m, tmp, x, y, j);

double * addFirstQuarter(double * m, double * quarter, unsigned int x, unsigned int y, unsigned int p ){

	unsigned int offsetM = 0, offsetQ = 0, offsetTMP =0;
	double* tmp = new double[x*y*4];

	for(int i = 0; i < y*2 ; i++){

		memcpy(tmp + offsetTMP, m + offsetM, sizeof(double) * x*2);
		offsetTMP += x*2;
		offsetM += sqrt(p);

	}

	offsetTMP = 0;

	for(int i = 0; i < y ; i++){

		memcpy(tmp + offsetTMP, quarter + offsetQ, sizeof(double) * x);
		offsetTMP += x*2;
		offsetQ += x;

	}

	//std::cout << "la" <<std::endl;
	//ecrit_bmp256("quarter.bmp", x*2, x*2, tmp);
	return tmp;
}

void amr2D_97(double* m, int p, int j) {

	double* final = new double[p] ;
	double* tmp = new double[p];
	memcpy(tmp, m, sizeof(double) * p);

	unsigned int x = 0, y = 0;
	x = sqrt(p);
	y = sqrt(p);

	for(unsigned int n=0; n < j; n++){

		analyse2D_97(tmp, x*y);

		writeFirstQuarter(final, tmp, x, y, static_cast<int>(sqrt(p)));

		tmp = extractFirstQuarter(tmp, x, y);

		x = static_cast<int>(x/2);
		y = static_cast<int>(y/2);

	}
	delete[] tmp;

	ecrit_bmp256("lena_AMR2D_3niv.bmp", sqrt(p), sqrt(p), final);
	memcpy(m, final, sizeof(double) * p);
	delete[] final;


}

void iamr2D_97(double* m, int p, unsigned int j){


	int size = p/(pow(2,(j-1)));

	double* tmp = new double[size];

	unsigned int x = 0, y = 0, offset = 0, offsetM = 0;
	x = sqrt(p)/pow(2, j-1);
	y = sqrt(p)/pow(2, j-1);

	//Copie du premier jeu de 4 quarts
	for(int i = 0; i < y; i++){

		memcpy(tmp + offset, m + offsetM, sizeof(double) * x);
		offset += x;
		offsetM += sqrt(p);

	}

	for(unsigned int n=0; n < j; n++){

		synthese2D_97(tmp, x*y);

		if(n <j-1) tmp = addFirstQuarter(m, tmp, x, y, p);

		x = static_cast<int>(x*2);
		y = static_cast<int>(y*2);


	}
	ecrit_bmp256("lena_IAMR2D_3niv.bmp", sqrt(p), sqrt(p), tmp);



	memcpy(m, tmp, sizeof(double) * p);
	delete[] tmp;


}

// quart 1 = superieur droit
// quart 2 = inférieur droit
// quart 3 = inférieur gauche
void statSousBande(double* m, unsigned int p,  unsigned int quart , unsigned int j, std::vector<double>* varList){


	int x = sqrt(p)/(pow(2,j));
	int y = sqrt(p)/(pow(2,j));

	unsigned int offsetTMP = 0, offsetM = 0;

	double* tmp = new double[x*y];

	switch(quart){

	case 0:{

		offsetM = 0;
				for(int i = 0; i < y; i++){

					memcpy(tmp + offsetTMP, m + offsetM, sizeof(double) * x);
					offsetTMP += x;
					offsetM += sqrt(p);

				}
				//ecrit_bmp256("quart.bmp", x, y, tmp);
				break;
	}

	case 1: {

		offsetM = x;
		for(int i = 0; i < y; i++){

			memcpy(tmp + offsetTMP, m + offsetM, sizeof(double) * x);
			offsetTMP += x;
			offsetM += sqrt(p);

		}
		//ecrit_bmp256("quart.bmp", x, y, tmp);
		break;
	}
	case 2: {


		offsetM = p/pow(2,j) + x;
		for(int i = 0; i < y; i++){

			memcpy(tmp + offsetTMP, m + offsetM, sizeof(double) * x);
			offsetTMP += x;
			offsetM += sqrt(p);

		}
		//ecrit_bmp256("quart2.bmp", x, y, tmp);

		break;
	}
	case 3: {

		offsetM = p/pow(2,j);
		for(int i = 0; i < y; i++){

			memcpy(tmp + offsetTMP, m + offsetM, sizeof(double) * x);
			offsetTMP += x;
			offsetM += sqrt(p);

		}
		//ecrit_bmp256("quart.bmp", x, y, tmp);

		break;
	}

	default: {
		break;
	}

	}//end switch

	//calcul esperance
	double espCarr =0, carrEsp = 0, var = 0;

	for(int i =0; i<y; i++){
		for(int j =0; j<x; j++){

			espCarr += tmp[i*y+j]*tmp[i*y+j];
			carrEsp	+= tmp[i*y+j];

			}
	}

	espCarr /= x*y;
	carrEsp /= x*y;
	double esp = carrEsp;
	carrEsp *=	carrEsp;

	var = espCarr - carrEsp;

	//Cherche min

	double min = 99999999;
	for(int i =0; i<y; i++){
			for(int j =0; j<x; j++){

				if(tmp[i*y+j] < min) min = tmp[i*y+j];

				}
		}

	//std::cout << "Var " << var << " min " << min << std::endl;
	FILE* file;
	file = fopen("statSousBande.txt", "a");
	fprintf (file, "Var sous bande  %d au niveau %d : %f \n", quart, j, var);
	fprintf (file, "Min sous bande  %d au niveau %d : %f \n", quart, j, min);
	fprintf (file, "Moy sous bande  %d au niveau %d : %f \n\n", quart, j, esp);
	fclose(file);
	varList->push_back(var);
	delete[] tmp;


}

void ajusterCoef(double* m, int p, int niveau){



	unsigned int x = sqrt(p)/pow(2, niveau);
	unsigned int y = sqrt(p)/pow(2, niveau);


	double* coefApprox = new double[x*y];
	unsigned int offset = 0, offsetM = 0;
	for(int i = 0; i< y; i++){

		memcpy(coefApprox + offset, m + offsetM, sizeof(double)*x );
		offset += x;
		offsetM += static_cast<int>(sqrt(p));

	}

	double max = maxSousBande(coefApprox, x*y);
	double coef = max/255;

	//Ajustement de la dynamique des coef d'approximation
	for(int i = 0; i < x*y; i++){

		coefApprox[i] /= coef;
	}

	//Ajustement des coef de detail
	for(int i = 0; i < p; i++){

		m[i] *= 8;
		m[i] += 10;
	}

	//Copie du coef d'approximation dans limage globale
	offset = 0, offsetM = 0;
	for(int i = 0; i < y; i++){

		memcpy(m + offsetM, coefApprox + offset , sizeof(double) * x );
		offset += x;
		offsetM += sqrt(p);

	}

	ecrit_bmp256("ajust.bmp", sqrt(p), sqrt(p), m);
	std::cout << "max " << max << std::endl;


}


void calculDebit(std::vector<double>* varList , int i, double b){

	double produit = 0;
	double result = 0;
	double sigma_i = varList->at(i);


	produit = puissance(varList->at(9) , 64*64);// sb 0 niv 3
	produit *= puissance(varList->at(8), 64*64);// sb 3 niv 3
	produit *= puissance(varList->at(7), 64*64);// sb 2 niv 3
	produit *= puissance(varList->at(6), 64*64);// sb 1 niv 3
	produit *= puissance(varList->at(5), 128*128);// sb 3 niv 2
	produit *= puissance(varList->at(4), 128*128);// sb 2 niv 2
	produit *= puissance(varList->at(3), 128*128);// sb 1 niv 2
	produit *= puissance(varList->at(2), 256*256);// sb 3 niv 1
	produit *= puissance(varList->at(1), 256*256);// sb 2 niv 1
	produit *= puissance(varList->at(0) , 256*256);// sb 1 niv 1


	result = sigma_i / produit;
	result = log(result)/log(2);
	result *= 0.5;

	result += b;

	std::cout << "result " << result << std::endl;
	FILE *file;
	file = fopen("debitSousBande.txt", "a");
	fprintf (file, "debit binaires %f :  \n", result);
	fclose(file);

}

double puissance(double sigma, double Nj){

	double result =0;

	double NJ_sur_N = (Nj / (512*512));
	result = pow(sigma, NJ_sur_N);
	//std::cout << "pow test" << result << std::endl;
//	result = pow(sigma, NJ_sur_N);

	return result;


}
